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Solitons confined in a channel are studied in the two-dimensional nonlinear Schrodinger equation. 
When a channel branches into two channels, a soliton is split into two solitons, if the initial kinetic 
energy exceeds a critical value. The branching point works as a pulse splitter. If it is below the 
critical value, the soliton is reflected. The critical kinetic energy for splitting is evaluated by a 
variational method. The variational method can be applied in the design of channel systems with 

f^S ■ reduced reflection. 
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The one-dimensional nonlinear Schrodinger equation is a typical soliton equation. Optical solitons in optical fibers 
are described by the nonlinear Schrodinger equation. [l[ The Bose- Einstein condensates (BECs) can be described by 
the Gross-Pitaevskii equation, which is equivalent to the nonlinear Schrodinger equation with a potential term.[2[ 
Bright and dark solitons have been observed in several experiments of the Bosc-Einstein condensates. [3J, |j] The one- 
dimensional nonlinear Schrodinger equation has been sufficiently studied, but solitons in two- or three-dimensional 
■ nonlinear Schrodinger equations are not so intensively studied. We studied solitons in guided channels in the two- 
\ dimensional nonlinear Schrodinger equation. Q In optical systems, refraction-index-guiding channels are used to con- 
p I . fine optical pulses. Matter-wave solitons in BECs are created in cigar-shaped traps, which work as guiding channels. 
Solitons can propagate along the guiding channel with an arbitrary velocity if the guiding channel is uniform and the 
norm of the soliton is below a critical value for the collapse. 

In this paper, we study the motion of solitons in inhomogeneous guiding channels. In particular, we study the 
motion of a soliton when a channel branches into two channels. Such a branching point is important as a pulse 
splitter in optical systems and a beam splitter for atomic waves. @ There has been no direct numerical simulation 
of the splitting of two-dimensional solitons in the branching waveguide up to now. We will perform a variational 
approximation using the Lagrangian to understand the numerical results. The approximation is a generalization 
^vq applied to the one-dimensional nonlinear Schrodinger equation under an external potential. [7| 
. The model equation is written as 

£v ^ = -iv 2 0-|^| 2 ^ + C/(a ; , 2 /)^ (1) 

where U(x,y) denotes the potential, which represents a channel. For a straight channel, U{x,y) = —Uq for — xq < 
x < xq and U {x, y) — for other regions. The width of the channel is denoted as 2xo and the depth of the potential 
is denoted as Uo- A stationary solution to eq. (1) can be numerically obtained from the time evolution of the 
Ginzburg-Landau-type equation with a feedback term: 
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a(N - N), (2) 



where /i is an additional variable corresponding to the chemical potential, N = J J \(f>\ 2 dxdy is the total norm, and 
Nq is the target value of the total norm. By a long-time evolution of eq. (2), a two-dimensional soliton with norm 
No is obtained as an attractor of the dynamical system. Figure 1(a) displays a 3D plot of <j> for Nq = 5, Uo = 5, and 
xo = 3. Figure 1(b) shows the profile \4>{x p ,y)\ in the section of x = x p = 0, and Fig. 1(c) shows \<p(x,y p )\ in the 
section of y = y p , where (x p ,y p ) is the peak position of the soliton solution. 

The Lagrangian corresponding to the nonlinear Schrodinger equation (1) is expressed as 

L=f y [(l/2){z(6'0/6't)^* - i(^V^)^>} " (1/2)|V</>| 2 + 1/2|^.| 4 - C/(x, y)|^| 2 ]dxdy. (3) 

If (j> is approximated as <j> — Asech(y/a) exp{— x 2 / (2b 2 )} exp(—ifit), the Lagrangian is calculated as 
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FIG. 1: (a) 3D plot of \tj>\ at xo = 3, Uo = 5, and No — 5. (b) Profile of \<f>\ in the section of x = 0. The dashed curve is 
an approximate solution obtained by the variational method, which well overlaps the numerical result, (c) Profile of \<j>\ in the 
section of y — y p , where (0, y p ) is the peak position of shown in (a). The dashed curve is the approximate solution 

obtained by the variational method. 



where N — J J \<j)\ 2 dxdy — 2A 2 y* r Kab is the total norm and erf(x) = 2/^pK f exp(—z 2 )dz is the error function. The 
variational principle dL e /da = dL e /db — yields coupled equations for a and b as 

a = 2\f2^b/N 1 1/2- N 2 /(24n) + 2U / v ^enp{-xl/b 2 )(-x a b) = 0. (5) 

The approximate values of a — 1.357,6 = 1.353, and A = 0.876 are obtained by the variational method for 
N = 5, Uq = 5, and xo = 3. The dashed curve in Fig. 1(b) denotes <f> — Asech{(y — £)/a} where £ = 20, and the 
dashed curve in Fig. 1(c) denotes 4> = Aexp{—x 2 /(2b 2 )}. The variational method is a good approximation for the 
two-dimensional soliton in the channel. 

Next, we study a problem in which the channel width changes in the y-direction as one of the simplest inhomogeneous 
systems. We assume that the channel width Xo(y) becomes narrow as Xo(y) = £01 for y < y\, Xo(y) = X01 + (xo2 — 
Xoi) x (y — (2/2 — 2/i) f° r 2/i < y < 2/2: and xo(y) = X02 for y > y^. The center line of the channel is fixed to x — 0. 
We have performed direct numerical simulation from the initial condition cj)(x,y) = 4>o(x, y) exp(ifcy), where <f>o(x,y) 
is the stationary solution shown in Fig. 1 for N — 5,xq =3, and Uq = 5, and k is the initial wave number in the 
y direction. If the channel width is constant, the two-dimensional soliton propagates with the velocity v y — k. The 
channel width changes smoothly, since yi — y\ is much larger than x§\ — £02 • The numerical simulation was performed 
by the pseudospectral method with 128 x 1024 Fourier modes. The system size was 20 x 160. Figure 2(a) shows the 
time evolution for k = 0.77, Xqi — 3, 202 = 2, 2/1 = 40, and y-2 = 100. The initial peak position of the soliton is (0, 20). 
The soliton propagates into the narrower channel, but the velocity becomes slower and the soliton amplitude becomes 
higher in the narrower channel. Figure 2(b) shows the time evolution at k = 0.74 for the same channel. The soliton 
cannot penetrate into the narrow channel and it is reflected. The threshold wave number k c is evaluated as k c = 0.76 
for the channel. We have investigated the threshold value k c for various values of X02 for xoi = 3, Uo = 5, y% — 40 and 
t/2 = 100. The rhombi in Fig. 2(c) show the relation of k c vs X02 by direct numerical simulation. That is, k c decreases 
with xq2- 

To understand this phenomenon, we performed a variational method by assuming <p(x,y,t) = Ascch{(y — 
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FIG. 2: (a) Time evolution of \<f)(x,y,t) \ in the section of x — at k = 0.77. (b) Time evolution of \<j>(x,y,t)\ in the section of 
x — fat k = 0.74. (c) Critical value k c for reflection. The dashed curve denotes the critical curve estimated by the variational 
method. 
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£(i))/a} exp{— x 2 /(2b 2 )} exp{ip(y — — i[it}. For the assumption, the effective Lagrangian is expressed as 



L eff = N{pC t (t)-p 2 /2 
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where f/ e // — l/(4^ 2 ) + l/(6a 2 ) — N/(6>/2nab) — Uoer{(xo/b). The parameters a and b are evaluated using the 
variational principle at each point y = £(t). Equation (7) is equivalent to Newton's equation of motion with the 
effective potential U e ff(y), which increases monotonically between y\ and 7/2 • The threshold value k c of the wave 
number k is evaluated from k 2 /2 = U e ff(x 02 ) — U e ff(x i). The dashed curve in Fig. 2(c) shows the threshold value 
evaluated by the variational method. The agreement with the direct numerical simulation is fairly good, but the 
theoretical curve is slightly larger than the numerically obtained values. 

The main problem in this study is a problem in which one channel branches into two channels. The potential 
U(x,y) is assumed to be U(x,y) = —Uq for — xqi < x < xqi when y < yi, U(x,y) = —Uq for — X12 < x < —X13 and 
£13 < x < X12, U(x,y) = —Ui for — x\ 3 < x < X13 when y\ < y < yt, and U(x,y) = —Uq for — £02 < x < —X03 
and x 03 < x < x 02 when y > y 2 , where x 12 = x i + {x 02 - x i) x (y - yi)/(y 2 - yi), x 13 = x Q3 x (y - yi)/{y 2 - yi), 
and U\ = Uq x (y 2 — y)/(y 2 — y\). U(x,y) = in the other region. We introduced the Ui component in the 
potential U(x,y) for the potential to change smoothly in the y-direction. Figure 3(a) shows the channel region with 
U(x,y) = —Uq for xqi — 3, £03 = 3, £02 = 6,j/i = 40, and y 2 — 100. The initial condition is the same as before: 
4>(x,y,t) ~ Ascch{(y — 20)/a} exp{— £ 2 /(26 2 )} exp{ik(y — 20)}. Figure 3(b) shows the time evolution of \<j)(0,y,t)\ 
in the section of x = at k — 0.4. The other parameters are Uo — 5 and N = 5. The soliton is reflected at the 
branching point. Figure 3(c) shows the time evolution of \4>(x,y p )\ in the y-section including the peak position of the 
modulus \cf>\ at k — 0.8. The soliton with a single peak is split into two solitons propagating along the two channels. 
We investigated the threshold value k c for the transition from reflection to splitting, and obtained k c ~ 0.74. Figure 
3(d) shows the trajectory at the peak position for k = 0.8 in the (x, y) space. The peak position is located at x = 
for y < 43.8 and it is split into two peaks for y > 43.8. The peak positions continuously separate away from x = 0. 

For the variational method applied to this splitting process, (j> is assumed to be ^ = (A/2)sech{(y — 
£(i))/a}[exp{— (x — ri) 2 /(2b 2 )} + cxp{ — (x + rj) 2 / (2b 2 )}] exp{ip(y — £(i)) — ifit}. The variational principle yields, 
for this ansatz, 
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FIG. 3: (a) Branching channel with 101 = 3, X03 ~ 3, X02 = 6,2/1 = 40, and j/2 = 100. (b) Time evolution of \<j)(x,y,t)\ at the 
section of x — for k = 0.4. (c) Time evolution of \<f>(x, y, t)\ at the a;-section passing through the peak position (x v ,y v ) for 
k — 0.8. (d) Trajectory of the peak position (x p (t) , y p (t)) at k = 0.8. The dashed curve is obtained by the variational method. 
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where cci = 0, x 2 = xqi for y < yi, x\ — £13, and x 2 = £12 for V\ < V < J/2, and x\ — X03, and x 2 = xq 2 for y > y 2 . 
In Fig. 3(d), the peak positions evaluated by the variational method are drawn with dashed curves. The trajectory is 
well approximated by the variational method. The threshold value obtained by the variational method is k c = 0.77, 
which is comparable to the numerical one, k c ~ 0.74, but the theoretical estimate is slightly larger than the numerical 
one, which is similar to the case of the tapering channel. 

The reflection of the pulse at the branching point might be undesirable for applications to pulse splitting. We can 
design a system including a branching point with reduced reflection using the variational method. The difference 
U e ff(y 2 ) — U e ff(yi) of the effective potential was the origin of the reflection. The reflection is expected to disappear 
or be reduced, if the channel system is designed, in which the difference \U e ff(y 2 ) — U e ff(yi)\ is minimized. We can 
control both the channel width xo 2 — X03 and the depth of potential Uo(y). We have designed the value of Uo(y) using 
an evolution equation: dUo(y)/d,T — —j3(U e ff(yx) — U e ff(y)) with j3 > as U e ff(y) approaches J7 e //(j/i) at each point 
y > y±. The U\ component is assumed to be the same as before: U\{y) = 5x (y 2 — y)/(y 2 — y\). Figure 4(a) shows 
the designed value of Uo(y) for xoi — 3, X02 = 2, y\ = 40, and y 2 — 100. Figure 4(b) shows the time evolution of the 
soliton with the initial wave number k = 0.2. Even for this small wave number, the reflection is suppressed and the 
soliton is split into two solitons, that penetrate into the two channels. 

In summary, we have studied the motion of two-dimensional solitons confined in guided channels. We have found 
that solitons with small wave number are reflected at the branching point of the channel. The variational method 
can predict the reflection and splitting phenomena fairly well. The method can be applied in the design of channel 
systems with reduced reflection. 
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FIG. 4: (a) Depth of potential designed for no reflection by the variational method, (b) Time evolution of \<f>(x, y, i)\ for 
the designed potential in the rr-section passing through the peak position at k — 0.2. (c) Time evolution of \(f>(x,y,t)\ for the 
designed potential at k — 0.1 in the a;-section passing through the peak position. 



